Optimized ensemble Monte Carlo simulations of dense Lennard-Jones fluids 
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We apply the recently developed adaptive ensemble optimization technique to simulate dense 
Lennard-Jones fluids and a particle-solvent model by broad-histogram Monte Carlo techniques. 
Equilibration of the simulated fluid is improved by sampling an optimized histogram in radial 
coordinates that shifts statistical weight towards the entropic barriers between the shells of the 
liquid. Interstitial states in the vicinity of these barriers are identified with unprecedented accuracy 
by sharp signatures in the quickly converging histogram and measurements of the local diffusivity. 
The radial distribution function and potential of mean force are calculated to high precision. 
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The free energy landscapes of complex systems are 
characterized by many local minima that are separated 
by entropic barriers. The simulation of such systems 
with conventional Monte Carlo [lj or molecular simula- 
tion methods is slowed down by long relaxation times 
due to the suppressed tunneling through these barriers. 
Extended ensemble simulations address this problem by 
broadening the range of phase space which is sampled in 
the respective reaction coordinate. Recently, an adap- 
tive algorithm |3( was introduced that explores entropic 
barriers by sampling a broad histogram in a reaction co- 
ordinate and iteratively optimizes the simulated statisti- 
cal ensemble defined in the reaction coordinate to speed 
up equilibration. The key idea of the approach is to 
measure the local diffusivity along the reaction coordi- 
nate, thereby identifying the bottlenecks in the simula- 
tions and then using this information to systematically 
shift statistical weights towards these bottlenecks in a 
feedback loop. The so optimized histogram converges to 
a characteristic form exhibiting peaks at the bottlenecks 
of the simulation, e.g. in the vicinity of the entropic bar- 
riers. The simulation of an optimized ensemble results in 
equilibration times which can be substantially lower com- 
pared to other extended ensemble simulations that aim 
at samplinga flat histogram in the respective reaction 
coordinate 0, 3 . Flat- histogram algorithms include the 
multicanonical method p|, simulated and parallel tem- 
pering broad histograms Q and transition matrix 
Monte Carlo § when combined with entropic sampling, 
as well as the adaptive algorithm of Wang and Landau 
and subsequent improvements |10| and variations [TH 
thereof. 

In this paper we demonstrate how the ensemble opti- 
mization method can be applied to a system that exhibits 
multiple barriers along a continuous reaction coordinate. 
Specifically, we simulate dense Lennard-Jones (LJ) fluids 
in two dimensions with a pair-wise potential of the form 
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where e is the interaction strength and a a length param- 
eter. In the fluid the interaction between two particles 



at distance R in the presence of many surrounding parti- 
cles can be described by a potential of mean force (PMF) 
that has the form 
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where g(R) is the radial pair distribution function. The 
PMF and radial pair distribution function show a sinu- 
soidal modulation revealing shell structures that form 
in the liquid as illustrated for various temperatures in 
Fig-E Between the shells multiple entropic barriers form 
which confine particles to the shells and thereby suppress 
thermal equilibration between the shells. Over the last 
decades Lennard-Jones fluids have been extensively in- 
vestigated by numerical simulations 0,0 and provide an 
ideal test system to study novel algorithms. Extended 
ensemble simulations of dense Lennard-Jones fluids have 
been performed in energy space using Wang-Landau sam- 
pling |l2l and in radial coordinates using an adaptive 
integration method (AIM) |14j . Here we determine the 




FIG. 1: Potential of mean force for a Lennard-Jones fluid at 
density p — 0.8 for varying temperatures. The inset shows the 
vapor-liquid equilibrium (solid line) of a Lennard-Jones fluid 
calculated by an extended ensemble simulation in Ref. |12|. 
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optimized ensemble in radial coordinates and its char- 
acteristic histogram for various temperatures and sub- 
sequently calculate the radial distribution function and 
potential of mean force. Approaching vapor-liquid co- 
existence we observe distinct features in the optimized 
histogram that allow us to identify interstitial states in 
the vicinity of the entropic barriers between the shells of 
the liquid. 

Ensemble optimization - In our Monte Carlo simula- 
tions the system does a Markovian random walk in con- 
figuration space which is described by particle positions 
in real space. This random walk is projected onto a ran- 
dom walk in two reaction coordinates, namely onto the 
energy E of a configuration which can be calculated by 
evaluating the pair-wise potential in Eq. as well as 
onto the radial distance R of two arbitrarily chosen par- 
ticles in the system. In energy space we simulate a canon- 
ical ensemble with weights w(E) = exp(—(3E), where we 
fix the inverse temperature (3. For the second reaction 
coordinate, the radial distance R between two particles, 
we define an additional broad-histogram ensemble with 
weights w(R) as explained below. These two ensembles 
define the acceptance probabilities for moves based on 
the Metropolis scheme 
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FIG. 2: (Color online) Radial local diffusivity D(R) for a 
dense Lennard- Jones fluid for varying temperatures. The free 
energy barriers between the shells of the liquid (see lower 
panel) lead to a suppression of the local diffusivity. At 
low temperatures additional peaks (arrows) emerge reveal- 
ing interstitial states between the shells of the liquid. The 
dashed lines indicate the interstitial barriers to which addi- 
tional weight is shifted for the optimized ensemble/histogram. 
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For sufficiently long simulation times we can measure 
the equilibrium histogram in the radial ensemble n w (R) 
which for any ensemble w(R) allows us to determine the 
radial distribution function g(R) — n w (R)/w(R). 

To speed up equilibration between the shells we opti- 
mize the statistical ensemble following the adaptive al- 
gorithm described in Ref. Q. Initially we sample a flat 
histogram n w (R) — w(R)g(R) by setting the statistical 
weights to a fixed value w(R) — l/g(R), where a rough 
estimate of the radial distribution function g(R) was ob- 
tained by a short run of the Wang-Landau algorithm 9]- 
During the simulation batch we measure round-trips of 
the random walk along the radial reaction coordinate by 
adding a label to the random walker that indicates which 
of the two extremal radial distances i? m i n or i? m ax it has 
visited most recently. We record two histograms, n w (R) 
and n+(i?), where the histogram n w (R) is incremented 
in every Monte Carlo step, while the histogram n+(R) 
is increased only if the label indicates that with respect 
to the two extrema, i? m ; n and -R max , the random walker 
has visited i? m i n most recently. This allows to define a 
current from i? m i„ to R 

max 



j = D(R)n w (R) 
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where D{R) is the local diffusivity at distance R and 
f{R) = n+(R)/n w (R) is the fraction of walkers which 
have most recently visited R m i n . Rearranging Eq. Q 



gives a simple means to measure the local diffusivity 



D(/,>I ~ ' "■■ AH, 7Ti 



(5) 



In Ref. |3( it was derived that in order to speed up equili- 
bration along the reaction coordinate the current j can be 
maximized by sampling an optimized histogram that is 
inversely proportional to the root of the local diffusivity 
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To sample such a histogram the weights ui( OJ,t ) (i?) of an 
improved statistical ensemble are obtained as 



w (°pt)(R) = W (R). 
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which arc then used to define the acceptance probabili- 
ties in Eq. J3J for a subsequent simulation batch. This 
feedback loop is iterated by simulating batches with re- 
peatedly refined statistical weights w(R) and increasing 
statistical accuracy. In our implementation we double 
the number of MC sweeps in subsequent batches, where 
one sweep corresponds to one attempted update per par- 
ticle. The weights w(R) are found to converge within 
four feedback iterations and a total of some 10 6 — I0 7 
MC sweeps. 

Simulations - We performed simulations for a 2D 
Lennard-Jones fluid with fixed density p = 0.8 at vary- 
ing temperatures. We used a relatively small system of 
TV = 144 particles confined to a square box of length 
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FIG. 3: (Color online) Optimized histograms for a dense 
Lennard- Jones fluid after feedback of the local diffusivity for 
varying temperatures. For the optimized ensemble peaks 
evolve at the free energy barriers between the shells of the 
liquid which proliferate at lower temperatures. Additional 
peaks emerge at low temperatures (dashed lines) revealing 
interstitial states (arrows) between the shells of the liquid. 




MC sweeps 



FIG. 4: Average statistical error (A In g(R))n of the com- 
puted radial distribution function for the optimized ensem- 
ble and flat-histogram sampling. Results are averaged over 
f6 independent Monte Carlo runs. The inset shows the de- 
viation of the radial distribution function obtained for the 
optimized ensemble from a reference distribution function 
S(R) = In g(R) — lng TC {(R) for 16 independent runs with 10 8 
MC sweeps. The reference distribution function was obtained 
from high statistics runs. 



13.42cr with periodic boundary conditions. We simulated 
the pair-wise Lennard-Jones potential given in Eq. 
without using a cut-off or shift. The continuous radial 
coordinate R was binned into units of size O.Olcr. In every 
update step a random displacement 5 < O.lcr in x- and y- 
directions was suggested for a single particle and accepted 
according to the acceptance probability in Eq. In 
order to calculate the radial distribution function g(R) 
we simulated a broad histogram in the radial coordi- 
nate covering a range [i? m i n , i? ma x] ~ [0.95er, 6.0<t]. For 
each batch the local diffusivity is determined by record- 
ing two histograms, n w (R) and n+(i?), and calculating 
the derivative df/dR by a linear regression which then 
allows to perform the feedback step in Eq. Q). 

In a first step we measure the local diffusivity D(R) 
for the projected random walk along the radial distance 
R. The striking feature of the measured local diffusivity 
shown in Fig. [21 is its strong variation along this reaction 
coordinate. Specifically, one observes for high temper- 
atures T > 1 that the local diffusivity has a maximum 
at the locations of the shells of the liquid (correspond- 
ing to the minima in the PMF) and is suppressed where 
the entropic barriers form between the shells. After feed- 
back of the local diffusivity the sampled histogram is no 
longer flat, but peaks emerge at the locations of the bar- 
riers between the shells, see Fig. [3] From the sampled 
histogram we have calculated the potential of mean force 
as $pmf(-R) = (lnw(R) — \nn w (R)) / (3 which is shown 
in Fig. ^ We find that the numerical accuracy of the 
estimated radial distribution function and PMF is mod- 
erately improved with an average speedup of 50 — 100% 
by sampling the optimized histogram instead of a flat 
histogram as shown in Fig. 



Interstitial states - As the temperature is further low- 
ered to T < 1 and the system is driven towards vapor- 
liquid coexistence, see Fig. Q] additional features emerge 
in both the measured local diffusivity and the optimized 
histogram, indicated by the arrows in Figs. [21 and [3J that 
point to the formation of additional structures. In the lo- 
cal diffusivity a sharp additional peak emerges between 
the first and second shell at radial distance R « 1.94ct 
as shown in Fig. [21 The close proximity of this peak 
to R = \[%\/2a reveals that an interstitial state becomes 
stable at low temperatures which is formed when the first 
shells around the two particles merge as illustrated in 
Fig. O b). Similarly, additional peaks between the higher 
shells point to interstitial states that are formed when 
the second/third shells around two particles merge. 

The distinct signature of the interstitial states in these 
plots unveils the potential use of the optimized ensemble 



a) separate shells b) interstitial state c) merged shells 




FIG. 5: Shell structures in a dense Lennard-Jones fluid, a) 
Two particles at large distance with separate shells and a total 
of 12 neighbors, b) An interstitial state forms as the shells 
around two particles merge with a characteristic distance of 
R = V3V2a « 1.94cr and a total of 10 neighbors, c) Pair 
of particles confined to the first shell with a characteristic 
distance of R — \/2a ~ 1.12cr and a total of 8 neighbors. 
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FIG. 6: (Color online) Optimized histograms for a particle- 
solvent model with varying interaction strength e. The solid 
lines indicate barriers between the shells of the liquid, the 
dashed lines refer to interstitial barriers, and the arrows indi- 
cate interstitial states. The second panel shows the average 
force, that is the derivative of the potential of mean force 
shown in the third panel. 

Monte Carlo method to identify interstitial states in the 
simulation of more complex systems such as fluids un- 
dergoing^ liquid- liquid phase transition [l5| , amorphous 
states pl| or biological systems. For instance, recent neu- 
tron diffraction ex peri ments |l7j | with high-density amor- 
phous (HDA) ice [HI revealed an interstitial molecular 
state which distinguishes HDA from low-density amor- 
phous (LDA) ice. Molecular dynamics simulations 0] 
confirmed this scenario by measuring the radial oxygen- 
oxygen distribution function g(R). While the intersti- 
tial state can be easily tracked in the pair distribution 
function for the high-density phase, the signature quickly 
faints as the pressure is reduced and a more sensitive 
measure such as the optimized histogram may help to 
identify the nature of the transition. For comparison we 
have indicated the location of the interstitial states also 
in the potential of mean force and its derivative, the av- 
erage force, in Fig. [H] In contrast to the optimized his- 



togram only a weak feature in the average force points 
to the interstitial state between the first and second shell 
which was also pointed out in Ref. [T3 |. 

Particle-solvent model - As another application of the 
optimized ensemble method we have studied a simple 
particle solvent model where we vary the strength e given 
in Eq. £Q) of an attractive interaction between a pair of 
particles and a solvent. Within the solvent we assume 
uniform coupling e = 1 . Varying the interaction strength 
e between the particles and the solvent is equivalent to a 
local renormalization of the temperature. As the inter- 
action strength is increased at a fixed temperature T = 1 
shells of the solvent form around the particles as illus- 
trated in Fig. a) . When the two particles approach 
each other the break-up of these shells is protected by a 
free energy barrier. Optimizing the radial ensemble sta- 
tistical weight is shifted towards this barrier and a peak 
in the histogram emerges at the location of the barrier 
as shown in Fig. El For increasing interaction strength e 
this peak at R ~ 2.14a further proliferates and becomes 
the dominant feature in the histogram. Due to the lo- 
cal character of the temperature renormalization upon 
increasing e the peaks at the locations of interstitial bar- 
riers between higher shells are only weakly amplified. 

Conclusions - In conclusion, by simulating a dense 
Lennard-Jones liquid we demonstrated that the opti- 
mized broad-histogram ensemble for a system with multi- 
ple entropic barriers can be efficiently determined by the 
recently developed adaptive ensemble optimization tech- 
nique [2j. The algorithm allows to systematically shift 
statistical weight towards the barriers between the shells 
of the liquid by measuring the local diffusivity along a 
continuous radial reaction coordinate. Distinct features 
in the measured local diffusivity and optimized histogram 
point to interstitial states and allow to significantly in- 
crease the sensitivity to detect such states compared to 
other approaches based on measurements of the radial 
distribution function or average force. The ensemble 
optimization method thereby provides a simple means 
to detect interstitial states that should be applicable to 
more complex systems. 
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